package com.yqg.education.utils;

public class test {

    // CGCS2000椭球参数
    private static final double CGCS2000_A = 6378137.0; // 长半轴
    private static final double CGCS2000_F = 1 / 298.257222101; // 扁率
    private static final double CGCS2000_E2 = 2 * CGCS2000_F - CGCS2000_F * CGCS2000_F;

    // GCJ-02加密参数（克拉索夫斯基椭球）
    private static final double GCJ02_A = 6378245.0;
    private static final double GCJ02_EE = 0.00669342162296594323;

    public static void main(String[] args) {
//        double x = 4364376.266; // 北坐标（米）
//        double y = 37637217.28; // 东坐标（含带号）
        double x = 4364196.733; // 北坐标（米）
        double y = 37637228.19; // 东坐标（含带号）

        // 步骤1：解析带号和中央子午线
        int zone = (int) (y / 1000000);
        double centralMeridian = zone * 3; // 3度分带
        double easting = y - zone * 1000000 - 500000; // 实际东坐标（米）

        // 步骤2：高斯反算为经纬度
        double[] wgs84 = gaussToWGS84(x, easting, centralMeridian);
        double lat = wgs84[0];
        double lon = wgs84[1];

        // 步骤3：转换为GCJ-02
        double[] gcj02 = wgs84ToGcj02(lat, lon);
        System.out.printf("GCJ-02坐标: 纬度=%.6f, 经度=%.6f", gcj02[0], gcj02[1]);
    }

    /**
     * 高斯投影反算（CGCS2000椭球）
     * @param northing 北坐标（米）
     * @param easting  东坐标（米，已扣除带号和500km）
     * @param l0       中央子午线经度（度）
     * @return [纬度, 经度]（十进制度）
     */
    private static double[] gaussToWGS84(double northing, double easting, double l0) {
        double a = CGCS2000_A;
        double e2 = CGCS2000_E2;

        // 迭代计算底点纬度
        double Bf = northing / (a * (1 - e2/4 - 3*e2*e2/64 - 5*e2*e2*e2/256));
        double lastBf;
        do {
            lastBf = Bf;
            double M = a * ((1 - e2/4 - 3*e2*e2/64 - 5*e2*e2*e2/256) * Bf
                    - (3*e2/8 + 3*e2*e2/32 + 45*e2*e2*e2/1024) * Math.sin(2*Bf)
                    + (15*e2*e2/256 + 45*e2*e2*e2/1024) * Math.sin(4*Bf)
                    - (35*e2*e2*e2/3072) * Math.sin(6*Bf));
            Bf += (northing - M) / (a * (1 - e2/4 - 3*e2*e2/64 - 5*e2*e2*e2/256));
        } while (Math.abs(Bf - lastBf) > 1e-10);

        // 计算卯酉圈曲率半径
        double Nf = a / Math.sqrt(1 - e2 * Math.sin(Bf) * Math.sin(Bf));
        double tf = Math.tan(Bf);
        double eta2 = e2 / (1 - e2) * Math.cos(Bf) * Math.cos(Bf);

        // 计算纬度B和经差l
        double yNorm = easting / Nf;
        double B = Bf - tf/(2*Nf*Nf) * easting * easting
                + tf/(24*Nf*Nf*Nf*Nf) * (5 + 3*tf*tf + eta2 - 9*eta2*tf*tf) * Math.pow(easting,4);
        double l = yNorm - yNorm*yNorm*yNorm/6 * (1 + 2*tf*tf + eta2)
                + yNorm*yNorm*yNorm*yNorm*yNorm/120 * (5 + 28*tf*tf + 24*tf*tf*tf*tf);

        return new double[] {
                Math.toDegrees(B),
                l0 + Math.toDegrees(l)
        };
    }

    /**
     * WGS-84/CGCS2000转GCJ-02
     */
    private static double[] wgs84ToGcj02(double lat, double lon) {
        if (isOutOfChina(lat, lon)) return new double[]{lat, lon};

        double[] delta = calculateDelta(lat, lon);
        return new double[]{lat + delta[0], lon + delta[1]};
    }

    private static double[] calculateDelta(double lat, double lon) {
        double dLat = transformLat(lon - 105.0, lat - 35.0);
        double dLon = transformLon(lon - 105.0, lat - 35.0);

        double radLat = Math.toRadians(lat);
        double magic = Math.sin(radLat);
        magic = 1 - GCJ02_EE * magic * magic;
        double sqrtMagic = Math.sqrt(magic);

        dLat = (dLat * 180.0) / ((GCJ02_A * (1 - GCJ02_EE)) / (magic * sqrtMagic) * Math.PI);
        dLon = (dLon * 180.0) / (GCJ02_A / sqrtMagic * Math.cos(radLat) * Math.PI);
        return new double[]{dLat, dLon};
    }

    private static double transformLat(double x, double y) {
        double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * Math.sqrt(Math.abs(x));
        ret += (20.0 * Math.sin(6.0 * x * Math.PI) + 20.0 * Math.sin(2.0 * x * Math.PI)) * 2.0 / 3.0;
        ret += (20.0 * Math.sin(y * Math.PI) + 40.0 * Math.sin(y / 3.0 * Math.PI)) * 2.0 / 3.0;
        ret += (160.0 * Math.sin(y / 12.0 * Math.PI) + 320 * Math.sin(y * Math.PI / 30.0)) * 2.0 / 3.0;
        return ret;
    }

    private static double transformLon(double x, double y) {
        double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * Math.sqrt(Math.abs(x));
        ret += (20.0 * Math.sin(6.0 * x * Math.PI) + 20.0 * Math.sin(2.0 * x * Math.PI)) * 2.0 / 3.0;
        ret += (20.0 * Math.sin(x * Math.PI) + 40.0 * Math.sin(x / 3.0 * Math.PI)) * 2.0 / 3.0;
        ret += (150.0 * Math.sin(x / 12.0 * Math.PI) + 300.0 * Math.sin(x / 30.0 * Math.PI)) * 2.0 / 3.0;
        return ret;
    }

    private static boolean isOutOfChina(double lat, double lon) {
        return lon < 72.004 || lon > 137.8347 || lat < 0.8293 || lat > 55.8271;
    }
}
